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I. INTRODUCTION 

What does dynamical systems theory contribute to our 
understanding of matter? To a large extent, the royal 
road to gain an understanding of fluids or solids has 
been statistical mechanics. Based on interaction poten¬ 
tials obtained from experiments and quantum mechani¬ 
cal simulations, sophisticated perturbation theories are 
capable of providing a quantitative description of the 
static structural properties of such fluids on the atom¬ 
istic scale. Computer simulations have provided guidance 
and invaluable insight to unravel the intricate local struc¬ 
ture and even the non-equilibrium dynamics of “simple” 
liquids including hydrodynamic flows and shock waves. 
At present, this program is being extended to ever more 
complex molecular fluids and/or to systems confined to 
particular geometries such as interfaces and pores, and 
to fluids out of thermal equilibrium. 

We have raised a related question, “What is liquid?” 
in 1998 in a similar context [T]. A fluid differs from a solid 
by the mobility of its particles, and this ability to flow 
is a collective phenomenon. The flow spreads rapidly in 
phase space, which constitutes the fundamental instabil¬ 
ity characteristic of a gas or liquid. About thirty years 
ago, dynamical systems theory provided new tools for 
the characterization of the chaotic evolution in a multi¬ 
dimensional phase space, which were readily applied to 
liquids shortly after mm- The main idea is to follow 
not only the evolution of a state in phase space but, si¬ 
multaneously, the evolution of various tiny perturbations 
applied to that state at an initial time and to measure 
the growth or decay of such perturbations. The study 
of the Lyapunov stability, or instability, of a trajectory 
with respect to small perturbations always present in real 
physical systems is hoped to provide new and alternative 
insight into the theoretical foundation and limitation of 
dynamical models for dense gases and liquids, of phase 
transitions involving rare events, and of the working of 
the Second Law of thermodynamics for stationary non¬ 
equilibrium systems. It is the purpose of this review to 
assess how far this hope has become true and how much 
our understanding of fluids has gained from the study of 
the Lyapunov instability in phase space. 


The structure of a simple fluid is essentially determined 
by the steep (e.g. oc repulsive part of the pair 

potential, which can be well approximated by a discon¬ 
tinuous hard core. In perturbation theories, this hard 
potential may be taken as a reference potential with the 
long-range attractive potential (oc —r“®) acting as a per¬ 
turbation [3]. Hard disk fluids in two dimensions, and 
hard sphere systems in three, are easy to simulate and 
are paradigms for simple fluids. The first molecular dy¬ 
namics simulations for hard sphere systems were carried 
out by Alder and Wainwright [S] in 1957, the first compu¬ 
tation of Lyapunov spectra for such systems by Dellago, 
Posch and Hoover [5] in 1996. 

There exist numerous extensions of the hard-sphere 
model to include rotational degrees of freedom and var¬ 
ious molecular geometries [7]. Arguably the simplest, 
which still preserves spherical symmetry, are spheres 
with a rough surface, so-called ’’rough hard spheres” 
[8]. In another approach, fused dumbbell diatomics are 
used to model linear molecules O [10] . Both models 
are used to study the energy exchange between trans¬ 
lational and rotational degrees of freedom and, hence, 
rotation-translation coupling for molecules with moder¬ 
ate anisotropy. Other approaches more suitable for larger 
molecular anisotropies involve spherocylinders smu], 
ellipsoids [131 |T3] and even inflexible hard lines [H] . To 
our knowledge, only for the first two schemes, namely for 
rough hard disks min] and for planar hard dumbbells 
[IKIHKIS], extensive studies of the Lyapunov spectra have 
been carried out. We shall come back to this work below. 

There are numerous papers for repulsive soft-potential 
systems and for Lennard-Jones systems from various au¬ 
thors, in which an analysis of the Lyapunov instability 
has been carried out. We shall make reference to some 
of this work in the following sections. 

The paper is organized as follows. In the next section 
we introduce global and local (time dependent) Lyapunov 
exponents and review some of their properties. Of par¬ 
ticular interest are the symmetry properties of the local 
exponents and of the corresponding perturbation vectors 
in tangent space. Both the familiar orthonormal Gram- 
Schmidt vectors and the covariant Lyapunov vectors are 
considered. In Sec. [m] we study planar hard disk sys- 
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terns over a wide range of densities, and compare them 
to analogous fluids interacting with a smooth repulsive 
potential. There we also demonstrate that the largest 
(absolute) Lyapunov exponents are generated by pertur¬ 
bations, which are strongly localized in physical space: 
only a small cluster of particles contributes at any in¬ 
stant of time. This localization persists in the thermo¬ 
dynamic limit. Sec. |IV| is devoted to another property of 
the perturbations associated with the small exponents, 
the so-called Lyapunov modes. In a certain sense, the 
Lyapunov modes are analogous to the Goldstone modes 
of fluctuating hydrodynamics (such as the familiar sound 
and heat modes). However, it is surprisingly difficult to 
establish a connection between these two different view¬ 
points. In Sec. |V] particle systems with translational and 
rotational dynamics are considered. Two simple planar 
model fluids are compared, namely gases of rough hard 
disks and of hard-dumbbell molecules. Close similarities, 
but also surprising differences, are found. Most surpris¬ 
ing is the fast disappearance of the Lyapunov modes and 
the breakdown of the symplectic symmetry of the local 
Gram-Schmidt exponents for the rough hard disk sys¬ 
tems, if the moment of inertia is increased from zero. In 
Sec. IVII we summarize some of the results for station¬ 
ary systems far from thermal equilibrium, for which the 
Lyapunov spectra have been valuable guides for our un¬ 
derstanding. For dynamically thermostatted systems in 
stationary non-equilibrium states, they provide a direct 
link with the Second Law of thermodynamics due to the 
presence of a multifractal phase-space distribution. We 
conclude with a few short remarks in Sec. m 


II. LYAPUNOV EXPONENTS AND 
PERTURBATION VECTORS 

Let T{t) = {p, q} denote the instantaneous state of 
a dynamical particle system with a phase space M of 
dimension D. Here, p and q stand for the momenta 
and positions of all the particles. The motion equations 
are usually written as a system of first-order differential 
equations, 

f = F(r), (1) 

where F is a (generally nonlinear) vector-valued function 
of dimension D. The formal solution of this equation is 
written as T{t) = 0*(r(O)), where the map </>* : F —F 
defines the flow in M, which maps a phase point F(0) at 
time zero to a point F(t) at time t. 

Next, we consider an arbitrary infinitesimal perturba¬ 
tion vector ST{t), which depends on the position F(t) in 
phase space and, hence, implicitly on time. It evolves 
according to the linearized equations of motion, 

,5r = J(F)(IF. (2) 

This equation may be formally solved according to 

6T{t) = Z?<(.‘|r(o) ^r(O), (3) 


where Dcj)* defines the flow in tangent space and is rep¬ 
resented by a real but generally non-symmetric D x D 
matrix . The dynamical (or Jacobian) matrix, 

OF 

^(r) ^ (4) 

determines, whether the perturbation vector i5F(t) has 
the tendency to grow or shrink at a particular point F(t) 
in phase space the system happens to be at time t. Ac¬ 
cordingly, the matrix element 

JFt(F),7(F)JF(F) (5) 

turns out to be positive or negative, respectively. Here, ^ 
means transposition. If, in addition, the perturbation is 
normalized, ||JF|| = 1, and points into particular direc¬ 
tions in tangent space to be specified below, this matrix 
element turns out to be a local rate for the growth or de¬ 
cay of I |i5F(t) 11 and will be referred to as a local Lyapunov 
exponent A(F) at the phase point F. 


A. Covariant Lyapunov vectors 


In 1968 Oseledec published his celebrated multiplica¬ 
tive ergodic theorem pPl - BH] . in which he proved that 
under very general assumptions (which apply to molecu¬ 
lar fluids) the tangent space decomposes into subspaces 
with dimension nij, 

TM{T) = (F) 0 (F) 0 ■ ■ • 0 (F), (6) 


for almost all F G M (with respect to the Lesbegue mea¬ 
sure), such that = D. These subspaces evolve 

according to 

Ino) , (7) 

and are said to be covariant, which means that they co¬ 
move - and in particular co-rotate - with the flow in 
tangent space. In general they are not orthogonal to 
each other. If v (F(0)) G (F(0)) is a vector in the 
subspace E^^'> (F(0)), it evolves according to 

7J<^‘|r(o) V (F(0)) = u (F(t)) G E^^^ (F(t)). (8) 


The numbers 


(±)y(i) 


lim -p-r 

i—>-±oo |t| 


\vim) 

vim) 


(9) 


for j G {1,2, ...,L} exist and are referred to as the 
(global) Lyapunov exponents. Here, the upper index (0) 
or (—) indicates, whether the trajectory is being followed 
forward or backward in time, respectively. If a subspace 
dimension rrij is larger than one, then the respective ex¬ 
ponent is called degenerate with multiplicity mj. If all 
exponents are repeated according to their multiplicity, 
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there are D exponents altogether, which are commonly 
ordered according to size, 


(+)Ai >•••>(+) Ai?, (10) 

(-)Ai >•••>(-) Ad, (11) 


where the subscripts are referred to as Lyapunov index. 
The vetors generating A^ according to 




1 , ||^'(r(i)) 

lim — In fr-r-— 
t^±oo |t| |ki^(r(o)) 


= 1,2,...,A 


( 12 ) 

are called covariant Lyapunov vectors. This notation, 
which treats the tangent space dynamics in terms of a set 
of vectors, is more convenient for algorithmic purposes, 
although the basic theoretical objects are the covariant 
subspaces ; j = 1, 2,..., L. 

Degenerate Lyapuov exponents appear, if there exist 
intrinsic continuous symmetries (such as invariance of 
the Lagrangian with respect to time and/or space trans¬ 
lation, giving rise to energy and/or momentum conserva¬ 
tion, respectively). For particle systems such symmetries 
almost always exist. Some consequences will be discussed 
below. 

The global exponents for systems evolving forward or 
backward in time are related according to 


=-(-^Ad+i-^; ee{l,2,...,D}. (13) 


The subspaces and, hence, the covariant Lyapunov 
vectors generally are not pairwise orthogonal. 

If in Eq. Q the time evolution is not followed from 
zero to infinity but only over a finite time interval r > 0, 
so-called covariant finite-time dependent Lyapunov expo¬ 
nents are obtained. 


(±)^pCOv = 1 In II D</±nr( 0 ) (r( 0 )) ||. ( 14 ) 


In the limit r —>■ 0 so-called local Lyapnnov exponents 
are generated, 


(±)ACOV(r) 


lim -p-r In II D(j)*\r (F) 
t^±o \t\ II ^ ^ ^ 


it 


(=^)t>^(r) j{T) (=^)t;^(r) 


(15) 


where ||v^(r)|| = 1 is required. The second equality has 
a structure as in Eq. (©, applied to the covariant vec¬ 
tors, and indicates that the local exponents are point 
functions, which only depend on F. The finite-time- 
dependent exponents may be viewed as time averages 
of local exponents over a stretch of trajectory during the 
finite time r, the global exponents as time averages over 
an infinitely-long trajectory. For the latter to become dy¬ 
namical invariants, one requires ergodicity, which we will 
always assume in the following for lack of other evidence. 


B. Orthogonal Lyapunov vectors 


Another definition of the Lyapunov exponents, also pi¬ 
oneered by Oseledec [201 - 123] , is via the real and symmet¬ 
rical matrices 


lim 

)-±oo 


D(j) 


tit 

lr(o) 


-D/'‘|r(o) 


2|t| 


(16) 


which exist with probability one (both forward and back¬ 
ward in time). Their (real) eigenvalues involve the global 
Lyapunov exponents, 

exp((=‘=^Ai) > • • • > exp((=‘=)AD), (17) 


where, as before, degenerate eigenvalues are repeated ac¬ 
cording to their multiplicities. For non-degenerate and 
degenerate eigenvalnes the > and = signs apply, re- 
spectively. The corresponding eigenspaces, Uf: ; j = 
1,... ,L, are pairwise orthogonal and provide two addi¬ 
tional decompositions of the tangent space at almost ev¬ 
ery point in phase space, one forward (-I-) and one back¬ 
ward (—) in time: 

TM(F) = t7i^^(F) 0 U^\r) 0 • ■ ■ 0 ui^\T). (18) 

In each case, the dimensions mj of the Oseledec sub¬ 
spaces t/± j = 1,... ,L, have a sum equal to the phase- 
space dimension D: These subspaces are 

not covariant. 

The classical algorithm for the computation of (global) 
Lyapunov exponents [MIES] carefully keeps track of all 
d-dimensional infinitesimal volume elements which 

(almost always) evolve according to 


(t) « (0) exp ( ^ Xet 


(19) 


This is algorithmically achieved with the help of a 
set of perturbation vectors ^^^g^(F); i = 1,...,D, 
which are either periodically re-orthonormalized with 
a Gram-Schmidt (GS) procedure (equivalent to a QR- 
decomposition) [27], or are continuously constrained to 
stay orthonormal with a triangular matrix of Lagrange 
multipliers HIIEHSI]- For historical reasons we refer to 
them as GS-vectors. They have been shown to be the 
spanning vectors for the Oseledec subspaces [32] . 
Accordingly, they are orthonormal but not covariant. 

Analogous to Eqs. (141 and (15), one may define finite- 

time-dependent GS Lyapunov exponents and 

local GS exponents (F), where the latter are 

again point functions, 

(±)ApS (F) = [(±V(F)]^ ,7(F) (±)gpF). (20) 


As before, the finite-time and global GS exponents are 
time averages of the respective local exponents over a 
finite respective infinite trajectory mM- 
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FIG. 1. (Color online) Covariant and GS vectors for a point 
O on the Henon attractor (light-blue). The black line is a 
finite-time approximation of the (global) stable manifold. 


From all GS vectors only which is associ¬ 

ated with the maximum global exponent, evolves freely 
without constraints which might affect its orientation in 
tangent space. Therefore it agrees with (^)ui(r) and is 
also covariant. Also the corresponding local exponents 
agree for i = 1, but generally not for other i. However, 
the local covariant exponents may be computed from the 
local GS exponents, and vice versa, if the angles between 
them are known [31) . 

As an illustration of the relation between covariant 
Lyapunov vectors and orthonormal GS vectors, we con¬ 
sider a simple two-dimensional example, the Henon map 

m, 

Xn+i = a - x^ + byn , 


C. Symmetry properties of global and local 
exponents 


For ergodic systems the global exponents are averages 
of the local exponents over the natural measure in phase 
space and, according to the multiplicative ergodic theo¬ 
rem for chaotic systems, they do not depend on the met¬ 
ric and the norm one applies to the tangent vectors. Also 
the choice of the coordinate system (Cartesian or polar, 
for example) does not matter. For practical reasons the 
Euclidian (or L2) norm will be used throughout in the 
following. It also does not matter, whether covariant 
or Gram-Schmidt vectors are used for the computation. 
The global exponents are truly dynamical invariants. 

This is not the case for the local exponents. They de¬ 
pend on the norm and on the coordinate system. And 
they generally depend on the set of basis vectors, covari¬ 
ant or GS, as was mentioned before. Furthermore, some 
symmetry properties of the equations of motion are re¬ 
flected quite differently by the two local representations. 


• Local Gram-Schmidt exponents: During the con¬ 
struction of the GS vectors, the changes of the dif¬ 
ferential volume elements following Eq. (19) 

are used to compute the local exponents. If the 
total phase volume is conserved as is the case for 
time-independent Hamiltonian systems, the follow¬ 
ing sum rule holds for almost all F: 


5]ApS(r) = o. 

i=i 


( 21 ) 


In this symplectic case we can even say more. For 
each positive local GS exponent there exists a neg¬ 
ative local GS exponent such that their pair sum 
vanishes [51] : 

(±)ApS(r) = -(±)AgS^_,(r), £ = 1,...,D. (22) 


2/n+l — j 

with a = 1.4 and b = 0.3. The light-blue line in Fig. [T] 
represents the Henon attractor, which is known to coin¬ 
cide with its unstable manifold. The black line is a finite¬ 
time approximation of its stable manifold. Let O denote 
a point on the attractor. The GS vectors at this point are 
indicated by QiiO) and 92 ( 0 ), where the former is also 
covariant and identical to vi{0) (parallel to the unsta¬ 
ble manifold). As required, the covariant vector V2(0) 
is tangent to the stable manifold (which was determined 
by a different method). For the computation of the co¬ 
variant vectors at O, it is necessary to follow and store 
the reference trajectory and the GS vectors sufficiently 
far into the future up to a point F in our example. In 
Fig. the GS vectors at this point are denoted by fl'i(E) 
and g 2 {F). Applying an algorithm by Ginelli et al. [55] . 
to which we shall return below, a backward iteration to 
the original phase point O yields the covariant vectors, 
V 2 { 0 ) in particular. 


As indicated, such a symplectic local pairing sym¬ 
metry is found both forward and backward in 
time. Non-symplectic systems do not display 
that symmetry. On the other hand, the re¬ 
orthonormalization process tampers with the ori¬ 
entation and rotation of the GS vectors and de¬ 
stroys all consequences of the time-reversal invari¬ 
ance property of the original motion equations. 


• Local covariant exponents: During their construc¬ 
tion [35j . only the norm of the covariant perturba¬ 
tion vectors needs to be periodically adjusted for 
practical reasons, but the angles between them re¬ 
main unchanged. This process effectively destroys 
all information concerning the d-dimensional vol¬ 
ume elements. Thus, no symmetries analogous to 
Eq. (22) exist. Instead, the re-normalized covariant 
vectors faithfully preserve the time-reversal symme¬ 
try of the equations of motion, which is reflected by 


(T) ACOV(r) = -(±)A«yi_,(r), £ = 1 ,..., D, (23) 
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regardless, whether the system is symplectic or not. 
This means that an expanding co-moving direction 
is converted into a collapsing co-moving direction 
by an application of the time-reversal operation. 

The set of all global Lyapunov exponents, ordered ac¬ 
cording to size, is referred to as the Lyapunov spectrum. 
For stationary Hamiltonian systems in thermal equilib¬ 
rium the (global) conjugate pairing rule holds. 


systems.. In all applications below appropriate reduced 
units will be used. 

To ease the notation, we shall in the following omit to 
indicate the forward-direction in time by (-I-), if there is 
no ambiguity. 

III. TWO-DIMENSIONAL HARD-DISK AND 
WCA FLUIDS IN EQUILIBRIUM 


(±)A. = 


D+l-f) 


which is a consequence of Eq. (22) and of the fact that the 


global exponents are the phase-space averages of these 
quantities. In such a case only the first half of the spec¬ 
trum containing the positive exponents needs to be com¬ 
puted. In the following we shall refer to this half as the 
positive branch of the spectrum. 

The subspaces spanned by the covariant (or GS) vec¬ 
tors associated with the strictly positive/negative global 
exponents are known as the unstable/stable manifolds. 
Both are covariant and are linked by the time-reversal 
transformation, which converts one into the other. 


D. Numerical considerations 

The computation of the Gram-Schmidt exponents is 
commonly carried out with the algorithm of Benettin 
et al. [24l [26] and Shimada et al. [25]. Based on 
this classical approach, a reasonably efficient algorithm 
for the computation of covariant exponents has been re¬ 
cently developed by Ginelli et al. |21] . Some com¬ 
putational details for this method may also be found in 
Refs. [1111311137]. An alternative approach is due to Wolfe 
and Samelson [33], which was subsequently applied to 
Hamiltonian systems with many degrees of freedom |39j . 

The considerations of the previous section are for time- 
continuous systems based on the differential equations 
Q. They may be readily extended to maps such as 
systems of hard spheres, for which a pre-collision state 
of two colliding particles is instantaneously mapped to 
an after-collision state. Between collisions the particles 
move smoothly. With the linearized collision map the 
time evolution of the perturbation vectors in tangent 
space may be constructed [5]. 

In numerical experiments of stationary systems, the 
initial orientation of the perturbation vectors is arbitrary. 
There exists a transient time, during which the pertur¬ 
bation vectors converge to their proper orientation on 
the attractor. All symmetry relations mentioned above 
refer to this well-converged state and exclude transient 
conditions [40]. 

For the computation of the full set of exponents, the 
reference trajectory and D perturbation vectors (each of 
dimension D) have to be followed in time, which requires 
D{D + 1) equations to be integrated. Present computer 
technology limits the number of particles to about 10^ 
for time-continuous systems, and to 10® for hard-body 


We are now in the position to apply this formalism to 
various models for simple fluids [41] . First we consider 
two moderately dense planar gases, namely a system of 
elastic smooth hard disks with a pair potential 




r < a, 
r > a, 


and a two-dimensional Weeks-Chandler-Anderson 
(WGA) gas interacting with a repulsive soft potential 

glia 


4>wca 



r < 2i/®cr, 
r > 2l/®cr. 


In Fig. [^the positive branches of their (global) Lyapunov 
spectra are compared. Both gases consist of = 400 
particles in a square box with side length L and periodic 
boundaries, and have the same density, p = N/L'^ = 0.4, 
and temperature, T = 1. Although, as expected, for low 
densities the Lyapunov spectra are rather insensitive to 
the interaction potential, the comparison of Fig. [^ for 



U2N 


FIG. 2. (Global) Lyapunov spectra of a planar hard-disk gas 
and of a planar WGA fluid with the same number of parti¬ 
cles, density and temperature. For details we refer to the main 
text. Only the positive branches of the spectra are shown. Al¬ 
though the spectra consist of discrete points for integer values 
of the index I, smooth lines and reduced indexes 1/2N are used 
for clarity in the main panel. In the inset a magnified view 
of the small-exponents regime is shown, for which I is not 
normalized. The figure is taken from Ref. m- 
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FIG. 3. Localization of the perturbation vector in physical 
space (the red square) for a planar gas of 102,400 soft repulsive 
disks |44| . The quantity g}, measuring the contribution of 
individual particles to the perturbation vector associated with 
the maximum GS exponent, is plotted in the vertical direction 
at the position of the particles. 



HAN 

FIG. 4. Localization spectra W for the complete set of Gram- 
Schmidt vectors (blue) and covariant vectors (red) for 198 
hard disks in a periodic box with an aspect ratio 2/11. The 
density p = 0.7 and the temperature T = 1. Reduced indices 
IjAN are used on the abscissa of the main panel. The in¬ 
set shows a magnihcation of the central part of the spectra 
dominated by Lyapunov modes. From Ref. m- 


moderately dense gases already reveals a rather strong 
sensitivity. In particular, the maximum exponent Ai is 
much lower for the WCA fluid, which means that de¬ 
terministic chaos is significantly reduced in systems with 
smooth interaction potentials. This difference becomes 
even more pronounced for larger densities, as will be 
shown below. 

The maximum (minimum) Lyapunov exponent de¬ 
notes the rate of fastest perturbation growth (decay) and 
is expected to be dominated by the fastest dynamical 
events such as a locally enhanced collision frequency. To 
prove this expectation one may ask how individual parti¬ 
cles contribute to the growth of the perturbations deter¬ 
mined by or . Writing the normalized perturbation 
vectors in terms of the position and momentum pertur¬ 
bations of all particles, f = 1,..., iV}, 

the quantity 

Mf'^(.Sq!'>)A(.Spf’)’ (24) 

is positive, bounded, and obeys the sum rule = 

1 for any i. It may be interpreted as a measure for the 
activity of particle i contributing to the perturbation in 
question. Equivalent relations apply for g^. If is 
plotted - vertical to the simulation plane - at the posi¬ 
tion of particle i, surfaces such as in Fig. are obtained. 
They are strongly peaked in a small domain of the physi¬ 
cal space indicating strong localization of = g^. It 
means that only a small fraction of all particles con¬ 
tributes to the perturbation growth at any instant of 
time. This property is very robust and persists in the 
thermodynamic limit in the sense that the fraction of 
particles, which contribute significantly to the formation 


of varies as a negative power of N [TFl mj 133 . For 
example, Fig. is obtained for a system of 102,400 (!) 
smooth disks interacting with a pair potential similar to 
the WCA potential [33]. 

This dynamical localization of the perturbation vectors 
associated with the large Lyapunov exponents is a very 
general phenomenon and has been also observed for one¬ 
dimensional models of space-time chaos [46] and Hamil¬ 
tonian lattices m- 

For £ > 1, the localization for differs from that of 
g^. This may be seen by using a localization measure due 
to Taniguchi and Morriss |3H1 EH] j 

^exp(5'); S{r{t)) = 


which is bounded according to 1/A^ ^ W < 1. The 
lower and upper bounds correspond to complete local¬ 
ization and delocalization, respectively. S is the Shan¬ 
non entropy for the ‘measure’ defined in Eq. (24), and 
(...) denotes a time average. The localization spectra 
for the two sets of perturbation vectors are shown in 
Fig. a for a rather dense system of hard disks in two 
dimensions. One may infer from the figure that localiza¬ 
tion is much stronger for the covariant vectors (red line) 
whose orientations in tangent space are only determined 
by the tangent flow and are not constrained by periodic 
re-orthogonalization steps as is the case for the Gram- 
Schmidt vectors (blue line). One further observes that 
the localization of the covariant vectors becomes less and 
less pronounced the more the regime of coherent Lya¬ 
punov modes, located in the center of the spectrum, is 
approached. 

Next we turn our attention to the maximum Lyapunov 
exponent Ai and to the Kolmogorov-Sinai entropy. Both 
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quantities are indicators of dynamical chaos. The KS- 
entropy is the rate with which information about an ini¬ 
tial state is generated, if a finite-precision measurement 
at some later time is retraced backward along the sta¬ 
ble manifold. According to Pesin’s theorem m it is 
equal to the sum of the positive Lyapunov exponents, 
= SAf>o determines the character¬ 

istic time for phase space mixing according to 

—mix _ ^ / h 

T — l/nKs- 

In Fig. we compare the isothermal density depen¬ 
dence of the maximum Lyapunov exponent Ai (top 
panel), of the smallest positive exponent X 2 N -3 (mid¬ 
dle panel), and of the Kolmogorov-Sinai (KS) entropy 
per particle Hks/N (bottom panel) for planar WCA and 
HD fluids. Both systems contain N = 375 particles in 
a periodic box with aspect ratio 0.6 and with a temper¬ 
ature T = 1. As expected, these quantities for the two 
fluids agree at small densities, but differ significantly for 
large p. 

For hard disks, van Zon and van Beijeren and de 
Wijn [55] used kinetic theory to obtain expressions for 
Ai and kxs/N to leading orders of p. They agree very 
well with computer simulations of Dellago et al. |56l [57] . 
The regime of larger densities, however, is only qualita¬ 
tively understood. and /N increase monoton- 

ically due to the increase of the collision frequency v. 
The (first-order) fluid-solid transition shows up as a step 
between the freezing point of the fluid ( pj^ = 0.88 ) 
and the melting point of the solid {p^^ = 0.91 |58|) (not 
shown in the figure). These steps disappear, if instead 
of the density the collision frequency is plotted on the 
abscissa. Af^^ and h^g /N diverge at the close-packed 
density due to the divergence of u. 

For the WCA fluid, both the maximum exponent and 
the Kolmogorov-Sinai entropy have a maximum as a 
function of density, and become very small when the den¬ 
sity of the freezing point is approached, which happens 
for pY^^ = 0.89 [55]. This behavior is not too surpris¬ 
ing in view of the fact that a harmonic solid is not even 
chaotic and Ai vanishes. The maximum of XY^^ occurs 
at a density of about 0.9pY^^ and confirms earlier re¬ 
sults for Lennard-Jones fluids m- At this density chaos 
is most pronounced possibly due to the onset of coop¬ 
erative dynamics characteristic of phase transitions. At 
about the same density, mixing becomes most effective 
as is demonstrated by the maximum of the KS-entropy. 
The latter is related to the mixing time in phase space 
according to = l/h^s [SUES]. 

It is interesting to compare the Lyapunov time = 
1/Ai, which is a measure for the system to “forget” its 
past, with the (independently-measured) time between 
collisions of a particle, Tc = l/i'. In Fig.j^such a compar¬ 
ison is shown for a three-dimensional system of 108 hard 
spheres in a cubic box with periodic boundary conditions 
[MUS]. Also included is the behavior of tks = N/kxs- 
For small densities, we have t\ <C Tc, and subsequent col¬ 
lisions are uncorrelated. This provides the basis for the 
validity of lowest-order kinetic theory (disregarding cor- 



FIG. 5. Isothermal density dependence of the maximum Lya¬ 
punov exponent, Ai (top), of the smallest positive exponent, 
X2N-3 (middle), and of the Kolmogorov-Sinai entropy per 
particle, kxs/N (bottom), for hard and soft-disk systems. 
From Ref. m- 


related collisions). For densities p > 0.4 the Lyapunov 
time is progressively larger than the collision time and 
higher-order corrections such as ring collisions become 
important. Also the lines for t\ and tks cross even be¬ 
fore, at a density 0.1. This is a consequence of the shape 
change of the Lyapunov spectrum with density |61] : the 
small positive exponents grow faster with p than the large 
exponents. It also follows that neither the Lyapunov time 
nor the mixing time determine the correlation decay of - 
say - the particle velocities. For lower densities, for which 





IV. LYAPUNOV MODES 



FIG. 6. (Color online) Comparison of the Lyapunov time t\ = 
1/Ai, the time tks = N/hxs, and the collision time Tc = \/v, 
as a function of the collision frequency v, for a system of 108 
hard spheres in a cubic box with periodic boundaries. The 
vertical lines indicate collision frequencies corresponding to 
the densities p = 0.1, 0.4, and 1.0. 


ring collisions are not important, the decay of the veloc¬ 
ity autocorrelation function z{t) is strictly dominated by 
the collision time. In Ref. [S3] we also demonstrate that 
for hard-particle systems the time t\ does not provide an 
upper bound for the time for which correlation functions 
may be reliably computed. 

For later reference, we show in Fig. the time- 
dependence of the local exponents for £ = 1 and £ = 
D = AN = 16 of a system consisting of four smooth hard 
disks. The symplectic symmetry as given by Eq. (221 is 
well obeyed. 



1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 
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FIG. 7. (Color online) Test of the symplectic symmetry of 
Eq. \22\ for a smooth hard-disk system with N = A. The 
phase space has 16 dimensions. The local GS exponents for 
£ = 1 and £ = 16 are plotted as a function of time t along a 
short stretch of trajectory. 


Already the first simulations of the Lyapunov spectra 
for hard disks revealed an interesting new step-like struc¬ 
ture of the small positive and negative exponents in the 
center of the spectrum [Bj, which is due to degenerate 
exponents. A similar structure was also found for hard 
dumbbell systems [Tj, to which we come back in Sec. 

The explanation for this behavior lies in the fact that 
the perturbation vectors for these exponents are charac¬ 
terized by coherent sinusoidal patterns spread out over 
the whole physical space (the simulation cell). We have 
referred to these collective patterns as Lyapunov modes. 
The modes are interpreted as a consequence of a sponta¬ 
neous breaking of continuous symmetries and, hence, are 
intimately connected with the zero modes spanning the 
central manifold [511 [B3|. They appear as sinusoidal mod¬ 
ulations of the zero modes in space - with wave number 
fc 0 - due to the spontaneous breaking of the contin¬ 
uous symmetries. For A: —>■ 0, the modes reduce to a 
linear superposition of the zero modes, and their Lya¬ 
punov exponent vanishes. The experimentally observed 
wave vectors depend on the size of the system and the 
nature of the boundaries (reflecting or periodic). 

Our discovery of the Lyapunov modes triggered a num¬ 
ber of theoretical approaches. Eckmann and Gat were 
the first to provide theoretical arguments for the ex¬ 
istence of the Lyapunov modes in transversal-invariant 
systems [64] . Their model did not have a dynamics in 
phase space but only an evolution matrix in tangent 
space, which was modeled as a product of independent 
random matrices. In another approach, McNamara and 
Mareschal isolated the six hydrodynamic fields related 
to the invariants of the binary particle collisions and the 
vanishing exponents, and used a generalized Enskog the¬ 
ory to derive hydrodynamic evolution equations for these 
fields. Their solutions are the Lyapunov modes [65] . In 
a more detailed extension of this work restricted to small 
densities, a generalized Boltzmann equation is used for 
the derivation |BB]. de Wijn and van Beijeren pointed out 
the analogy to the Goldstone mechanism of construct¬ 
ing (infinitesimal) excitations of the zero modes and the 
Goldstone modes of fluctuating hydrodynamics [BTlISBj. 
They used this analogy to derive the modes within the 
framework of kinetic theory [69] . With a few exceptions, 
a rather good agreement with the simulation results was 
achieved, at least for low densities. Finally, Taniguchi, 
Dettmann, and Morriss approached the problem from the 
point of view of periodic orbit theory m and master 
equations m- 

The modes were observed for hard-ball systems in one, 
two, and three dimensions miiia EH Ea na US], for 
planar hard dumbbells [18] . and also for one- and two- 
dimensional soft particles miEaiTg. 

A formal classification for smooth hard-disk systems 
has been given by Eckmann et al. m- If the AN com- 
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ponents of a tangent vector are arranged according to 


= {5qx,5qy\5px,5py), 


the six orthonormal zero modes, which span the central 
manifold and are the generators of the continuous sym¬ 
metry transformations, are given by iSlITSl 


ei 

62 

63 

64 

65 

66 


{Px,Py] 0,0), 

^( 1 , 0 ; 0 , 0 ), 
^( 0 . 1 ; 0 , 0 ), 

( 0 > 0 ; Px,Py), 
^( 0 , 0 ; 1 , 0 ), 
^( 0 , 0 ; 0 , 1 ). 


ei corresponds to a shift of the time origin, 64 to a change 
of energy, 62 and 63 to a uniform translation in the x and 
y directions, respectively, and 65 and eg to a shift of the 
total momentum in the x and y directions, respectively. 
The six vanishing Lyapunov exponents associated with 
these modes have indices 2N — 2 < i < 2N -|- 3 in the 
center of the spectrum. Since, for small-enough k, the 
perturbation components of a particle obey 5p = CJq 
for the unstable perturbations (A > 0), and (5q = —C6p 
for the stable perturbations (A < 0), where C > 0 is a 
number, one may restrict the classification to the 5q part 
of the modes and, hence, to the basis vectors 61 , 62,63 
m- Now the modes with a wave vector k may be clas¬ 
sified as follows: 

1. Transverse modes (T) are divergence-free vec¬ 
tor fields consisting of a superposition of sinusoidal 
modulations of 62 and 63 . 

2. Longitudinal modes (L) are irrotational vector 
fields consisting of a superposition of sinusoidal 
modulations of 62 and 63 . 


For the same k, the values of the Lyapunov expo¬ 
nents for the longitudinal and momentum modes coin¬ 
cide. These modes actually belong to a combined sub¬ 
space LP{n) = L(n) © P{n). The dimensions of the 
subspaces T{n) and LP{n) determine the multiplicity of 
the exponents associated with the T and LP modes. For 
a periodic boundary, say in x direction, the multiplicity 
is 2 for the T, and 4 for the LP modes. For details we 
refer to Ref. m- 

The transverse modes are stationary in space and time, 
but the L and P modes are not [73[7i[^. In the fol¬ 
lowing, we only consider the case of periodic boundaries. 
Any tangent vector from an LP subspace observed in a 
simulation is a combination of a pure L and a pure P 
mode. The dynamics of such an LP pair may be rep¬ 
resented as a rotation in a two-dimensional space with 
a constant angular frequency proportional to the wave 
number of the mode. 


= vkn, 

where the L mode is continuously transformed into the P 
mode and vice versa. Since the P mode is modulated with 
the velocity field of the particles, the experimentally ob¬ 
served mode pattern becomes periodically blurred when 
its character is predominantly P. Since during each half 
of a period, during which all spacial sine and cosine func¬ 
tions change sign, the mode is offset in the direction of the 
wave vector k, which gives it the appearance of a travel¬ 
ing wave with an averaged phase velocity v. If reflecting 
boundaries are used, standing waves are obtained. The 
phase velocity v is about one third of the sound velocity 
m, but otherwise seems to be unrelated to it. 

The basis vectors spanning the subspaces of the T and 
LP modes may be reconstructed from the experimental 
data with a least square method IS71I73]. As an example, 
we consider the L(1,0) mode of a system consisting of 198 
hard disks at a density 0.7 in a rectangular periodic box 
with an aspect ratio 2/11. The LP(1,0) subspace includes 
the tangent vectors for 388 < .^ < 391 with four identi¬ 
cal Lyapunov exponents. The mutually orthogonal ba¬ 
sis vectors spanning the corresponding four-dimensional 
subspace are viewed as vector fields in position space and 
are given by m 


3 . Momentum modes (P) are vector fields consist¬ 
ing of sinusoidal modulations of 64 . Due to the 
random orientations of the particle velocities, a P 
mode is not easily recognized as fundamental mode 
in a simulation. However, it may be numerically 
transformed into an easily recognizable periodic 
shape [mils!, as will be shown in Fig. below. 

The subspaces spanned by these modes are denoted by 
r(n), L(n), and P(n), respectively, where the wave 
vector for a periodic box with sides L^^Ly is k = 
(2'Knx/Lx,2'Kny/Ly), and n = (nx,riy). and Uy are 
integers. 



where, for simplicity, we have omitted the normalization. 
If the Gram-Schmidt vectors are used for the reconstruc¬ 
tion, the components corresponding to the cosine are 
shown for L(1,0) in Fig.[^ for P(1,0) in Fig.[^ In the lat¬ 
ter case, the mode structure is visible due to the division 
by Px and Py as indicated. This also explains the larger 
scattering of the points. The components proportional 
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FIG. 8. (Color online) Reconstructed L(1,0) mode compo¬ 
nents proportional to cos{2'Kqx/L^). The system consists of 
198 smooth hard disks in a periodic box. For details we refer 
to the main text. 

to the sine are analogous, but are not shown. A related 
analysis may be carried out for the covariant vectors. 

Recently, Morriss and coworkers CIHZil considered in 
detail the tangent-space dynamics of the zero modes and 
of the mode-forming perturbations over a time r, during 
which many collisions occur. In addition, they enforced 
the mutual orthogonality of the subspaces of the modes 
and of their conjugate pairs (for which the Lyapunov ex¬ 
ponents only differ by the sign), with the central manifold 
by invoking an (inverse) Gram-Schmidt procedure, which 
starts with the zero modes and works outward towards 
modes with larger exponents. With this procedure they 
were able to construct the modes in the limit of large r 
and to obtain (approximate) values for the Lyapunov ex¬ 
ponents for the first few Lyapunov modes (counted away 
from the null space) [78] . The agreement with simulation 
results is of the order of a few percent for the T modes, 
and slightly worse for the LP modes. 

For the modes with a wave length large compared 
to the inter-particle spacing, a continuum limit for the 
perturbation vectors leads to a set of partial differen¬ 
tial equations for the perturbation fields, whose solutions 


FIG. 9. (Golor online) Reconstrncted P(1,0) mode compo¬ 
nents proportional to cos{2TTqx/Lx) for the same LP snbspace 
as in Fig. 

give analytical expressions for the modes in accordance 
with the boundary conditions applied. This procedure 
works for the stationary transverse modes, as well as for 
the time-dependent LP modes. For the latter, the par¬ 
tial differential equations for the continuous perturbation 
fields assume the form of a wave equation with traveling 
wave solutions m For quasi-one-dimensional systems 
(with narrow boxes such that particles may not pass each 
other), the phase velocity is found to depend on the par¬ 
ticle density via the mean particle separation and the 
mean time between collisions (of a particle). For fully 
two-dimensional models of hard disks, the wave veloc¬ 
ity becomes proportional to the collision frequency of a 
particle, and inversely proportional to the mean squared 
distance of a particle to all its neighbors in the first co¬ 
ordination shell [771 izni- This is an interesting result, 
since it provides a long-looked-for connection between a 
mode property - the wave velocity - with other density- 
dependent microscopic quantities of a fluid. 
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V. ROTOTRANSLATION 
A. Rough hard disks 

Up to now we were concerned with simple fluids al¬ 
lowing only translational particle dynamics. As a next 
step toward more realistic models, we consider fluids, for 
which the particles may also rotate and exchange energy 
between translational and rotational degrees of freedom. 
In three dimensions, the first model of that kind, rough 
hard spheres, was already suggested by Bryan in 1894 
|80j . Arguably, it constitutes the simplest model of a 
molecular fluid. It was later treated by Pidduck m and, 
in particular, by Chapman and Cowling [5], who derived 
the collision map in phase space for two colliding spheres 
with maximum possible roughness. The latter requires 
that the relative surface velocity at the point of contact of 
the collision partners is reversed, leaving their combined 
(translational plus rotational) energy, and their combined 
linear and angular momenta invariant [7]. The thermo¬ 
dynamics and the dynamical properties of that models 
were extensively studied by O’Dell and Berne [SS]. They 
also considered generalizations to partial roughness in- 
between smooth and maximally-rough spheres |83H85j . 
Here we consider the simplest two-dimensional version of 
that model, maximally-rough hard disks in a planar box 
with periodic boundaries. Explicit expressions for the 
collision maps in phase space and in tangent space may 
be obtained from our previous work m mi [HS]. The 
coupling between translation and rotation is controlled 
by a single dimensionless parameter 

4/ 

K = - 

where I is the principal moment of inertia for a rotation 
around an axis through the center, and m and a denote 
the mass and diameter of a disk, respectvely. k may take 
values between zero and one: 0, if all the mass is concen¬ 
trated in the center, 1/2 for a uniform mass distribution, 
and 1, if all the mass is located on the perimeter of the 
disk. 

The rotation requires to add angular momentum Ji 
to the other independent variables qi,Pi of a particle 
i. The Ji show up in the collision map [811I61I86], and 
their perturbations 5Ji affect the tangent space dynamics 
[iniiHS]. If one is only interested in the global Lyapunov 
exponents, this is all what is needed, since it constitutes 
an autonomous system. We refer to it as the J-version 
of the rough-disk model. Its phase-space dimension D = 
5A. Since D may be an odd number, it is obvious that 
one cannot ask questions about the symplectic nature of 
the model. For this it is necessary to include also the 
disk orientations 0^ conjugate to the angular momenta 
in the list of independent variables. The angles do not 
show up in the collision map of the reference trajectory, 
but their perturbations, <50^, affect the collision map in 
tangent space. Now the phase space has 6N dimensions. 
We refer to this representation as the J0-version. 



FIG. 10. (Color online) Full Lyapunov spectra in J0- 
representation for four rough disks (ac = 0.5) respective 
smooth disks (ac = 0) with boundary conditions as indi¬ 
cated. It is interesting to note that the equipartition theo¬ 
rem does not hold for such small systems. For k = 0.5, the 
translational and rotational kinetic energies add according to 
3.59 -I- 2.41 = 6. For the simulation with ac = 0, a total trans¬ 
lational kinetic energy equal to two was used. 


Before entering the discussion of the Lyapunov spec¬ 
tra for large systems, let us clarify the number of van¬ 
ishing exponents first. As an illuminating case, we plot 
in Fig. [To] some spectra for the J0-version of a rough 
disk system with only 4 particles and 24 exponents. For 
positive K, when translation and rotation are intimately 
coupled, the number of vanishing exponents is deter¬ 
mined by the intrinsic symmetries of space and time- 
translation invariance and by the number of particles. 
For periodic boundaries in x and y directions (full red 
dots in Fig. 101, space-translation invariance contributes 
four, time-translation invariance another two vanishing 
exponents. Furthermore, the invariance of the collision 
map with respect to an arbitrary rotation of a parti¬ 
cle adds another N zero exponents, four in our exam¬ 
ple. Altogether, this gives 10 vanishing exponents, as 
is demonstrated in Fig. by the red points. For re¬ 
flecting boundaries in x and y directions (small blue dots 
in Fig. 101 there is no space-translation invariance, and 
the total number of vanishing exponents is limited to 
six. If K vanishes, the rotational energy also vanishes, 
and all variables connected with rotation cease to be 
chaotic and contribute another 2N zeroes to the spec¬ 
trum. Together with the intrinsic-symmetry contribu¬ 
tions, this amounts to 2N -|- 6 = 14 vanishing exponents 
for periodic boundaries (open purple circles in Fig. 101, 
respective 2N -|- 2 = 10 for reflecting boundaries (not 
shown). Note that we have plotted the full spectra in¬ 
cluding the positive and negative branches in this figure. 
For the J-version, and for an odd number of particles, the 
discussion is slightly more involved for which we refer to 
Ref. [I6]. 

Next we turn to the discussion of large systems, for 
which deviations from equipartition are negligible. We 
use reduced units for which m, a and the kinetic transla- 
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FIG. 11. (Color online) Lyapunov spectra (positive branches 
only) for a system of 88 rough hard disks for various coupling 
parameters k as indicated by the labels [16]. The density 
p = 0.7. The periodic box has an aspect ratio of 2/11. 



FIG. 12. (Color online) Enlargement of the regime of small 
Lyapunov exponents for some spectra of Fig. [m Also neg¬ 
ative exponents are displayed such that conjugate pairs, Ar 
and Asiv+i-r, share the same index I on the abscissa. 


tional (and rotational) temperatures are unity. Lyapunov 
exponents are given in units of \JK/Nma'^. The system 
we first consider in the following consists of TV = 88 rough 
hard disks with a density p = 0.7 in a periodic box with 
an aspect ratio A = 2/11 [T6|. For comparison, spec¬ 
tra for 400-particle systems are given in Ref. (SSj. Since 
we are here primarily interested in the global exponents, 
the simulations are carried out with the J-version of the 
model. 

In Fig. E] we show the positive branches of the Lya¬ 
punov spectra for a few selected values of k [TH|. Al¬ 
though the exponents are only defined for integer £, 
smooth lines are drawn for clarity, and a reduced in¬ 
dex is used on the abscissa. There are bN = 440 ex¬ 
ponents, of which the first half constitute the positive 
branch shown in the figure. An enlargement of the small 


exponents near the center of the spectrum is provided 
in Fig. 12 for selected values of k. There we have also 
included the negative branches in such a way to empha¬ 
size conjugate pairing. For k = 0, which corresponds 
to freely rotating smooth disks without roughness, the 
Lyapunov spectrum is identical to that of pure smooth 
disks without any rotation, if = 88 vanishing expo¬ 
nents are added to the spectrum. The simulation box 
is long enough (due to the small aspect ratio) for Lya¬ 
punov modes to develop along the x direction. This may 
be verified by the step structure of the spectrum due 
to the small degenerate exponents. If k is increased, the 
Lyapunov modes quickly disappear, and hardly any trace 
of them is left for k = 0.1. For small positive n the spec¬ 
trum is decomposed into a rotation-dominated regime 
for 2N < i < 3N and a translation-dominated regime 
for £ < 2N and i > 3N. It is, perhaps, surprising that, 
with the exception of the maximum exponent Ai, already 
a small admixture of rotational motion significantly re¬ 
duces the translation-dominated exponents, whereas the 
rotation-dominated exponents only gradually increase. 

To study this further, we plot in the Figs. and 
the isothermal K-dependence of the maximum exponent, 
Ai, and of the Kolmogorov-Sinai entropy per particle, 
^ks/N, respectively, for various fluid densities as indi¬ 
cated by the labels. These data are for a system con¬ 
sisting of 400 rough disks |HS|. Ai and, hence, dynamical 
chaos tends to decrease with k for low densities, and al¬ 
ways increases for large densities. At the same time, the 
KS-entropy always decreases with k. This means that 
mixing becomes less effective and the time for phase- 
space mixing goes up, the more the rotational degrees 
of freedom interfere with the translational dynamics. 

Before leaving the rough hard-disk case, a slightly dis¬ 
turbing observation for that system should be mentioned. 
In Fig. we have demonstrated that the local GS Lya¬ 
punov exponents of the smooth hard-disk model display 
symplectic symmetry. Unexpectedly, for the rough hard- 



K 


FIG. 13. (Color online) Dependence of the maximum Lya¬ 
punov exponent on the coupling parameter k for the rough- 
disk fluid at various densities. From Ref. |86| . 
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FIG. 14. (Color online) Dependence of the Kolmogorov-Sinai 
entropy per particle on the coupling parameter k for the 
rough-disk fluid at various densities. From Ref. m 



t 


FIG. 15. (Color online) Demonstration of the lack of symplec- 
tic symmetry (as formulated in Eq. ( |22[ )) for a rough hard¬ 
disk system with = 4 particles. The simulation is carried 
out with the J0-representation of the model, which includes 
the disk orientations, and for which the phase space has 24 
dimensions. The coupling parameter k = 0.5. The local GS- 
exponents for 1 = 1 and I = 2A are plotted as a function of 
time t for a short stretch of trajectory. 


B. Hard dumbbells 

A slightly more realistic model for linear molecules are 
“smooth” hard dumbbells, whose geometry is shown in 
Fig. As with hard disks, the dynamics is charac¬ 
terized by hard encounters with impulsive forces per¬ 
pendicular to the surface. Between collisions, the par¬ 
ticles translate and rotate freely, which makes the sim¬ 
ulation rather efficient [SlIinilHT]. In the following we 
restrict to the planar version of this model, for which 
the ’’molecule” consists of two fused disks as shown in 
Fig. The state space is spanned by the center-of- 
mass (CM) coordinates qi, the translational momenta Pi, 
the orientation angles 0^, and the angular momenta Ji, 
i = 1,2,...,A^, and has &N dimensions. The Lyapunov 
spectra for this model were hrst computed by Milanovic 
et al. pQ HSl [ini [EH] • Related work for rigid diatomic 
molecules interacting with soft repulsive site-site poten¬ 
tials was carried out by Borzsak et al. |89j and Kum et 
al. |90) . 

Here we restrict the discussion to homogeneous dumb¬ 
bells, for which the molecular mass m is uniformly dis¬ 
tributed over the union of the two disks. The moment of 
inertia for rotation around the center of mass becomes 

< [’’’ + 2 arctan(d/u>)] 

4 2(iw-I-cr2[7 r-I-2 arctan(d/u>)] 

and 


I{d > a) 


m 



where w = V— (P is the molecular waist. I{d) 
monotonously increases with d. For d —>■ 0 it converges 
to that of a single disk with mass m and diameter cr, and 
with a homogeneous mass density. Results for other mass 
distributions may be found in Ref. |18] . d plays a sim¬ 
ilar role as the coupling parameter k for the rough-disk 
model. 

The results are summarized in Fig. where the Lya¬ 
punov spectra, positive branches only, for a system of 


disks this is not the case. This is shown in Fig. 15 where 
the maximum and minimum local exponents do not show 
the expected symmetry, and neither do the other conju¬ 
gate pairs which are not included in the hgure. If the 
collision map for two colliding particles in tangent space 
is written in matrix form, dF' = AldF, where dF and 5F' 
are the perturbation vectors immediately before and after 
a collision, respectively, the matrix M does not obey the 
symplectic condition = M. Here, S is the anti¬ 

symmetric symplectic matrix, and ^ means transposition. 
For a more extensive discussion we refer to Refs. [Mini. 



FIG. 16. Geometry of a hard dumbbell diatomic. All quan¬ 
tities for this model are given in terms of reduced quantities, 
for which the disk diameter a, the molecular mass m, and the 
total kinetic energy per molecule, K/N, are equal to unity. 
The molecular anisotropy is defined by d/a. 
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FIG. 17. Anisotropy dependence of the Lyapunov spectra 
for a system of 64 planar hard dumbbells in a square box 
with periodic boundaries, d is the molecular anisotropy (in 
reduced units). Only the positive branches of the spectra with 
3N = 192 exponents are shown. The number density is 0.5. 
The Lyapunov index is denoted by 1. (From Ref. [T]). 


64 dumbbells in a periodic square box at an intermedi¬ 
ate gas density of 0.5 are shown for various molecular 
anisotropies d. The system is too small for Lyapunov 
modes to be clearly visible (although a single step due to 
a transverse mode may be observed for d > 0.01 even in 
this case). Most conspicuous, however, is the widening 
gap in the spectra between the indices I = 2N — 3 = 125 
and I = 126 for d < dc ~ 0.063. It separates the 
translation-dominated exponents (1 < ^ < 125) from the 
rotation-dominated exponents (126 < I < 189). The gap 
disappears for d> de¬ 
ll was observed in Ref. [T] that for d < dc the per¬ 
turbation vectors associated with the translation- and 
rotation-dominated exponents predominantly point into 
the subspace belonging to center-of-mass translation and 
molecular rotation, respectively, and very rarely rotate 
into a direction belonging to the other subspace. For 
anisotropies d > dc, however, one finds that the offset 
vectors for all exponents spend comparable times in both 
subspaces, which is taken as an indication of strong trans¬ 
lation - rotation coupling, dc is found to increase with 
the density. The anisotropy dependence of some selected 
Lyapunov exponents is shown in Fig. |18| The horizon¬ 
tal lines for / = 1 and I = 2N — 3 = 125 indicate the 
values of the maximum and smallest positive exponents 
for a system of 64 hard disks with the same density, to 
which the respective exponents of the dumbbell system 
converge with d —>■ 0. The smooth hard-disk data were 
taken from Ref. [5]. 

It was shown in the Refs. HUH] that two “phase tran¬ 
sitions” exist for the hard dumbbells with an anisotropy 
of d = 0.25, a first at a number density ni = 0.75 from a 
fluid to a rotationally-disordered solid, and a second, at 
712 ~ 0.775, from the orientationally-disordered solid to 
a crystal with long-range orientational order. Both tran¬ 



FIG. 18. Anisotropy dependence of selected Lyapunov expo¬ 
nents, labelled by their indices I, for the 64-dumbbell system 
summarized in Fig. |17| The horizontal lines indicate the val¬ 
ues for corresponding exponents of a smooth hard-disk gas. 
From Ref. [Tj. 


sitions were observed by computing orientational corre¬ 
lation functions. The first transition at ni makes itself 
felt by large fluctuations of the Lyapunov exponents dur¬ 
ing the simulation and a slow convergence [55], but does 
not show up in the density dependence of the Lyapunov 
exponents. The second transition at n 2 does give rise to 
steps in the density-dependence curves of both Ai and 
of Ai 89 (which is the smallest positive exponent). These 
steps are even more noticeable when viewed as a function 
of the collision frequency instead of the density, and the 
step is significantly larger for Aigg than for Ai. This in¬ 
dicates that the collective long-wavelength perturbations 
are much stronger affected by the locking of dumbbell 
orientations than the large exponents, which - due to 
localization - measure only localized dynamical events. 
Here, dynamical systems theory offers a new road to an 
understanding of collective phenomena and phase tran¬ 
sitions, which has not been fully exploited yet. 

Unlike the rough hard disks, the hard dumbbell sys¬ 
tems generate well-developed Lyapunov modes for their 
small exponents. Instructive examples for N = 400 ho¬ 
mogeneous hard dumbbells with an anisotropy d = 0.25 
are given in the Refs. [H] US]. The classification of 
the modes in terms of transversal (T) and longitudinal- 
momentum (LP) modes is completely analogous to that 
for hard disks m, although some modifications due to 
the presence of rotational degrees of freedom have not 
yet been worked out in all detail. In Fig. 19 we show, as 
an example, the smallest positive exponents of a system 
of = 150 dumbbells with a density n = 0.5 in an elon¬ 
gated (along the x axis) periodic box with an aspect ratio 
A = 1/32 [55] . Only modes with wave vectors parallel to 
X develop, which facilitates the analysis. The spectrum 
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FIG. 19. Magnification of the lower part of the Lyapunov 
spectrum for a system consisting of = 150 homogeneous 
dumbbells with a molecular anisotropy d = 0.25. The density 
n = 0.5, and the aspect ratio of the box is 1/32. The Lya¬ 
punov exponents are given in units of \/JjiJW)]~{rna^, where 
K is the total (translational and rotational) kinetic energy, m 
is the dumbbell mass, and cr is the diameter of the two rigidly 
connected disks forming a dumbbell. The complete spectrum 
contains 900 exponents. 
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FIG. 21. As in the previous Fig. |20| but for the longitudinal- 
momentum mode belonging to the perturbation vector (5r443. 


VI. STATIONARY PARTICLE SYSTEMS OUT 
OF EQUILIBRIUM 


displays the typical step structure due to the mode de¬ 
generacies with multiplicity two, for the T modes, and 
four, for the LP modes. For this system, typical instan¬ 
taneous mode patterns are shown in the Figs. [20| and [21] 
the former for a T mode with I = 446, the latter for an 
LP mode with I = 443. In both cases, the wavelength is 
equal to the box length in x direction. Please note that 
the system appears strongly contracted along the x axis. 



FIG. 20. Transverse mode pattern for the perturbation vector 
5r446 for an instantaneous configuration of the 150-dumbbells 
system with the Lyapunov spectrum shown in Fig. |19| The 
circles denote the center-of-mass positions of the dumbbells, 
and the arrows indicate the position perturbations (5®, 5y) 
(top) and momentum perturbations {Spx,Spy) (bottom) of 
the particles. The system is strongly contracted along the 
x-axis. 


One of the most interesting - and historically first - 
studies of Lyapunov spectra for many-particle systems 
is for stationary non-equilibrium states laElllII. If a 
system is driven away from equilibrium by an external 
or thermal force, the irreversibly generated heat needs to 
be removed and transferred to a heat bath. Otherwise 
the system would heat up and would never reach a steady 
state. Such a scheme may assume the following equations 
of motion, 

q, = Pi/rui + Q({q}, {p})A(t) 

p, = -d<^>/d<i, + T’({q}, {p})X(t) - SiCP», (25) 

where QX and VX represent the driving, and —CPi the 
thermostat or heat bath acting on a particle i selected 
with the switch Si S {0,1}. C fluctuates and may be pos¬ 
itive or negative, whether kinetic energy has to be ex¬ 
tracted from - or added to - the system. Averaged over 
a long trajectory, {() needs to be positive. C({q}i{p}) 
may be a Lagrange multiplier, which minimizes the ther¬ 
mostatting force according to Gauss’ principle of least 
constraint and keeps either the kinetic energy or the to¬ 
tal internal energy a constant of the motion. Or it may 
be a single independent variable in an extended phase 
space such as for the Nose-Hoover thermostat. There are 
excellent monographs describing these schemes [SH US] . 
They are rather flexible and allow for any number of 
heat baths. Although it is not possible to construct such 
thermostats in the laboratory, they allow very efficient 
non equilibrium molecular dynamics (NEMD) simula¬ 
tions and are believed to provide an accurate description 
of the nonlinear transport involved |23j . 

Simple considerations lead to the following thermody¬ 
namic relations for the stationary non-equilibrium state 
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generated in this way after initial transients have de¬ 
cayed: 



Here, d is the dimension of the physical space, and d Si 
becomes the number of thermostatted degrees of free¬ 
dom. D is the phase space dimension, and dQ/dt is the 
rate with which heat is transferred to the bath and gives 
rise to a positive rate of entropy production, (dSirr/dt). 
T is the kinetic temperature enforced by the thermostat. 
These relations show that an infinitesimal phase volume 
shrinks with an averaged rate, which is determined by the 
rate of heat transfer to the bath, which, in turn, is reg¬ 
ulated by the thermostat variable C- As a consequence, 
there exists a multifractal attractor in phase space, to 
which all trajectories converge for t —>■ oo, and on which 
the natural measure resides. For Axiom A flows or even 
non-uniform hyperbolic systems, these states have been 
identihed with the celebrated Sinai-Ruelle-Bowen (SRB) 
states |2S]. 

This result is quite general for dynamically thermostat¬ 
ted systems and has been numerically confirmed for many 
transport processes in many-body systems. The com¬ 
putation of the Lyapunov spectrum is essential for the 
understanding of this mechanism mEiiii] and, still, is 
the only practical way to compute the information di¬ 
mension of the underlying attractor. Such a system is 
special in the sense that its stationary measure is sin¬ 
gular and resides on a fractal attractor with a vanishing 
phase space volume. The ensuing Gibbs entropy is di¬ 
verging to minus infinity, which explains the paradoxical 
situation that heat is continuously transferred to the bath 
giving rise to a positive rate of irreversible entropy pro¬ 
duction Sirr ■ The fact that there is a preferred direction 
for the heat flow is in accordance with the Second Law of 
thermodynamics, which these systems clearly obey. The 
divergence of the Gibbs entropy also indicates that the 
thermodynamic entropy, namely the derivative of the in¬ 
ternal energy with respect to the entropy, is not defined 
for such stationary non-equilibrium states. 

We demonstrate such a fractal attractor for a one¬ 
dimensional Frenkel-Kontorova conductivity model of a 
charged particle in a sinusoidal potential, which is sub¬ 
jected to a constant applied field X and a Nose-Hoover 
thermostat [sniini]. The equations of motion are 

q = p-, p = - sm{q) + X - Cp; C=p^-1. 

The phase space is three-dimensional. In Fig. three 
Poincare maps for mutually orthogonal planes dehning 
the first octant of the phase space are shown. The sin¬ 
gularity spectrum for this model has been computed in 
Ref. [5D], and the multifractal nature of the stationary 
phase-space measure has been established. 





FIG. 22. (Color online) One-dimensional Frenkel-Kontorova 
conductivity model described in the main text. Three 
Poincare maps for the planes defining the first octant are 
shown. The system is in a stationary non-equilibrium state. 
Reprinted from Ref. [94] 


The sum of all Lyapunov exponents becomes negative 
for stationary non-equilibrium states. To demonstrate 
this effect, we show in Fig. [^the Lyapunov spectra of 
a system of 36 particles subjected to planar shear flow 
for various shear rates i as indicated by the labels. Ho¬ 
mogeneous SLLOD mechanics |92] in combination with a 
Gaussian energy control has been used, which keeps the 
internal energy of the fluid a constant of the motion. The 
figure emphasizes the conjugate pairing symmetry in the 
non equilibrium stationary state: Since all particles are 
homogeneously affected by the same thermostat, the pair 



FIG. 23. (Color online) Lyapunov spectra of a system of 36 
hard disks subjected to iso-energetic planar shear flow with 
shear rates e as indicated by the points. For details we refer 
to Ref. |6|. 
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sum of conjugate exponents becomes negative and is the 
same for all pairs. With the Kaplan-Yorke formula, the 
information dimension is found to be smaller than the 
dimension of the phase space [3], a clear indication of a 
multi-fractal phase-space distribution. 


VII. FINAL REMARKS 

Here we have described some of our studies of simple 
fluids, the results so far obtained, and the new questions 
they raise. Most of our simulations are for planar systems 
due to the numerical effort required for the computation 
of Lyapunov spectra, but the results are expected to carry 
over to three-dimensional systems. 

What have these studies taught us about fluids? For 
equilibrium systems, the results about the density depen¬ 
dence of the maximum Lyapunov exponents and of the 
mixing time in Sec. |III| clearly provide deeper and quanti¬ 
tative insight into the foundation and limitations of more 
familiar models for dense fluids. Still, the behavior of the 
exponents near phase transitions has not been satisfacto¬ 
rily exploited yet and should be the topic for further re¬ 
search. It is interesting to note that time-dependent local 
Lyapunov exponents have been used to identify atypical 
trajectories with very high or very low chaoticity. Such 
trajectories, although rare and difficult to locate in phase 
space, may be of great physical significance. For exam¬ 
ple, in phase traditions, or in chemical reactions, they 
may lead over unstable saddle points connecting a quasi¬ 
stable state, or chemical species, with another. Tailleur 
and Kurchan invented a very promising algorithm, Lya¬ 
punov weighted dynamics^ for that purpose [95] , in which 
a swarm of phase points evolves according to the natural 
dynamics. Periodically, trajectories are killed, or cloned, 
with a rate determined by their local exponents Ai. As a 
consequence, trajectories are weighted according to their 
chaoticity. A slightly modified version of this algorithm, 
Lyapunov weighted path sampling, has been suggested by 
Geiger and Dellago (Mj, and has been successfully ap¬ 
plied to the isomerizing dynamics of a double-well dimer 
in a solvent. 

Since the turn of the century, Lyapunov modes have 
raised a lot of interest. But it has been surprisingly dif¬ 
ficult to connect them to other hydrodynamic properties 


[TIITTlIZSj. Most recently, however, some progress has 
been made. But this problem will still occupy researchers 
for some time. 


Another interesting and promising field is the study 
of systems with qualitatively different degrees of free¬ 
dom, such as translation and rotation. Our studies of 
the Lyapunov instability of rough disks and of linear 
dumbbell molecules has opened another road to investi¬ 
gate translation-rotation coupling. An extension of these 
studies to vibrating molecules and to more complex sys¬ 
tems with internal rotation, such as rotational isomer 
dynamics in butane m, are promising topics for future 
research. 

Arguably, the most important impact Lyapunov spec¬ 
tra have achieved is in non-equilibrium statistical me¬ 
chanics. In combination with dynamical thermostats, it 
is possible to generate stationary states far from equi¬ 
librium and to link the rate of entropy production with 
the time-averaged friction coefficients and the logarith¬ 
mic rate of phase-volume contraction. See Eq.(26). The 
existence of a multi-fractal phase-space distribution with 
an information dimension smaller than the phase space 
dimension provides a geometrical interpretation of the 
Second Law of thermodynamics. The importance of this 
result is reflected in the voluminous literature on this 
subject. 

There are many other applications of local Lyapunov 
exponents and tangent vectors in meteorology, in the ge¬ 
ological science, and in other fields. For more details we 
refer to Ref. [M], which is a recent collection of review 
articles and applications. 
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